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1.0  INTRODUCTION 

Terrain  roughness  has  a  direct  impact  on  ground  vehicle  performance.  Movement  over  terrain  causes 
excitation  of  the  vehicle,  through  its  suspension  system,  based  on  forward  speed.  Its  effects  can  require 
reductions  in  speed  due  to  vibration  imparted  to  the  occupants  or  concern  about  structural  integrity. 
Performance  requirements  that  incorporate  terrain  roughness,  as  a  parameter,  ensure  vehicle  operation 
at  expected  levels  since  the  vehicle  will  be  designed  to  handle  excitations  for  that  environment. 
Accomplishment  of  this  requires  a  metric  that  measures  terrain  roughness  which  allows  vehicle  system 
designers  and  testers  the  ability  to  use  it  to  correlate  roughness  levels  with  associated  predicted  and 
verified  performance.  As  can  be  seen  from  references  [1]  to  [1 1]  this  has  been  a  field  of  interest  for  some 
time. 

The  scope  here  is  limited  to  the  metric  of  root-mean-square  (RMS)  as  a  measure  of  terrain  roughness — 
recognizing  its  limitations  (assumes  no  “favored  or  predominate  frequencies”)  and  accepting  its  historical 
significance  to  this  point  as  the  Army’s  measure  of  terrain  roughness. 


2.0  SPACIAL  DOMAIN 

The  description  and  examples  of  the  metric  are  given  using  MATLAB  [13]  as  the  computational  platform. 

2.1  Root-Mean-Square  (RMS)  Terrain  Roughness  Metric 

The  statistical  measure  of  standard  deviation  (equivalent  to  the  square  root  of  the  variance)  “shows  how 
much  variation  or  dispersion  from  the  average  exists”  (http://en.wikipedia.org/wiki/Standard  deviation). 
Another  term  for  this  is  root-mean-square  (RMS)  when  the  mathematical  operations  used  to  calculate  it 
are  considered.  The  assumption  that  an  average  exists  leads  to  the  requirement  that  the  metric  must  be 
used  on  stationary  data  to  make  sense.  As  Figure  1  shows,  the  RMS  roughness  metric  for  a  non¬ 
stationary  (non-detrended)  terrain  would  result  in  a  much  larger  value  than  for  its  stationary  (detrended) 
counterpart  and  not  provide  the  intended  correlation  between  vehicle  performance — the  RMS  “OF”  the 
trend  is  much  larger  than  the  RMS  of  the  profile  “FROM”  the  trend. 
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Figure  A-9  EFFECT  OF  DETRENDING  ON  DATA  SET  4 
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Figure  1  -  Stationary  Data  Example  (see  [12]) 


This  highlights  the  importance  of  detrending  before  calculating  terrain  roughness.  Leaving  the  issue  of 
which  detrending  method  is  best,  further  details  of  detrending  with  a  double-sided  exponential  filter  will  be 
discussed  in  the  formal  definition  of  the  RMS  roughness  metric. 

Energy  conservation  indicates  that  there  is  a  relationship  between  the  spacial  and  frequency 
representation  of  a  signal  (see  Parseval’s  theorem).  The  RMS  roughness  metric  of  a  terrain  profile 
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(spacial  domain)  is  an  equivalent  measure  of  the  area  under  the  Power  Spectral  Density  (PSD)  of  the 
same  profile  (frequency  domain).  From  [12]  the  terrain  PSD  can  be  defined  as 

PSDd(A )  =  CA~N 

(D 

where  ”N  is  approximately  2  for  both  natural  and  man-made  surfaces.  Man-made  surfaces  can  be 
artificially  constructed  to  give  any  value  of  N”.  With  the  definition  in  equation  (2) 


v 


(2) 


A 

wavelength  (distance/cycle) 

V 

velocity  (distance/time) 

f 

frequency  (cycle/time) 

the  dependence  of  the  excitation  of  the  vehicle  on  forward  velocity  is  shown  in  equation  (3). 

PSDd(f)  =  vCf~2 


(3). 


From  the  above  discussion,  the  formal  definition  of  the  RMS  terrain  roughness  metric  is  given  in  Figure  2. 


RMS  = 


ZJLoWi) 


N 


Fd&d  =  F(.Xi)  -  F(Xi) 


rna\ 

[  T7  f  'v  _L  -vi  /nr  _L  17  C  v  _  -ki  /nr  1  /o  l  3  ] 


F(x) 

elevation  at  point  x 

FdO) 

detrended  elevation  point 

F(x) 

trend  elevation  point 

n 

step  number 

a 

measurement  interval 

X 

weight  constant 
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M 


limit  of  summation 


Figure  2  -  RMS  Terrain  Roughness  Metric  Definition  (see  [12]) 

Application  of  the  metric  to  a  data  sequence  is  simple.  Making  the  data  stationary  (detrending)  to  obtain 
meaningful  results  is  the  difficulty — 

The  complexity  arises  when  the  principles  are  reduced  to  practice,  for  then  the  realities  of 
implementing  logic,  storing  information,  and  transferring  data  have  to  be  addressed. 

---Robert  F.  Stengel,  Optimal  Control  and  Estimation  (ISBN  0-486-68200-5) 


2.2  Double-sided  Exponential  Filter 

Based  on  work  done  in  [12],  a  double-sided  exponential  filter  has  historically  been  used  to  detrend  terrain 
elevation  profile  data  before  determining  its  RMS  value— the  necessity  of  which  has  been  discussed 
above.  In  electrical  engineering  terms  this  is  equivalent  to  high-pass  filtering  with  the  difference  that  for 
digital  data  double-sided  averaging  (both  past  and  future)  can  be  performed  to  remove  phase  shift  that 
may  result  from  one-sided  filtering. 

From  Figure  2  values  for  the  limit  of  summation,  measurement  interval,  and  weight  constant  must  be 
determined  to  apply  the  detrending  filter  as  shown  in  Figure  3. 


figure  A- 2  EXPONENTIAL  DETRENDING  Of  DATA  POINTS 


Figure  3  -  Double-sided  Exponential  Filter  (see  [12]) 

2.2.1  Measurement  Interval 

Conclusions  based  on  analysis  in  [12]  indicate  that  the  smallest  wavelength  that  affects  vehicles  from 
natural  terrains  is  about  2  feet.  With  consideration  for  the  Nyquist-Shannon  sampling  theorem  a 
maximum  of  1  foot  intervals  can  be  chosen. 


2.2.2  Limit  of  Summation 

Historically,  sixty  feet  (60  feet)  has  been  considered  the  longest  wavelength  that  affects  vehicle 
responses.  Equation  (4)  provides  an  example  assuming  the  lowest  sprung  mass  natural  frequency  is 
0.75  Hz  with  a  speed  of  30  MPH  (44  feet/sec). 
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v  44 

Aotai  =  j  =  ^5  =  58.6  =>  60  feet 


(4) 


This  is  the  complete  length  of  the  double-sided  exponential  filter — halving  the  value  for  one-side  (30  feet). 


2.2.3  Weight  Constant 

As  shown  in  equation  (5)  the  practical  numerical  limit  of  the  exponential  filter’s  extent  is  considered  met 
when  the  exponent  reaches  -3.  Including  terms  with  larger  exponents  will  contribute  little  since  their 
weighting  coefficients  will  be  quite  small. 


e  (I)  =  e-3  =  0.0497 

Table  1  summaries  the  filter  parameters  for  several  cases. 

_ Table  1  -  Exponential  Filter  Parameters 


exp  length  (1  side) 

30  ft 

360  in 

360  in 

a 

1  ft 

3  in 

6  in 

M  (-  len°th ) 

30 

120 

60 

/  Ma\ 

;(=ir) 

10ft 

120  in 

120  in 

(5) 


The  extent  is  3X  with  X  being  fixed  for  a  given  exponential  length  (1-side) — although  the  units  may  change 
(e.g.  feet  to  inches). 

2.3  Example  Usage  -  RMS  Terrain  Roughness  Metric  Function 

The  RMS  definition  in  Figure  2  is  implemented  in  Appendix  A.1  using  MATLAB  [13].  The  function  “rms”  is 
used  on  digital  elevation  versus  displacement  samples  in  the  spacial  domain.  The  assumption  for 
displacement  is  equally  spaced  data  samples. 

Example  -  calling  rms  function _ 

%  y  hold  digital  elevation  versus  displacement  data 

x  =  y(:,2)*12;  %  convert  feet  to  inches  or  consistent  units  for  other  input  values 
%  e- (N*A/lambda) 

N  =  [0:1:120]  %  1-side  =  30  ft  =  360  in  @  3  in  samples  =  120  (360/3) 

A  =  3;  %  3  in  samples 

lambda  =  10*12;  %  10  ft  *  12  =  120  in 

%  -N*A/lambda  =  -120*3/120  =  -3  @  last  sample 

[ fdet, rms , fmean, xmean] =rms (x, N, A, lambda) ;  %  units  are  inches 


The  output  of  this  function  will  be  used  to  compare  to  frequency  domain  calculations. 

3.0  FREQUENCY  DOMAIN  MATHEMATICS 

Historically  the  RMS  Terrain  Roughness  Metric  has  been  calculated  in  the  spacial  domain.  Digital  Signal 
Processing  (DSP)  generally  views  filters  from  a  frequency  domain  perspective.  Addressing  the  frequency 
domain  calculation  of  the  metric  in  Section  2.1  and  the  filter  defined  in  Section  2.2  should  result  in  wider 
application  and  understanding. 
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The  link  between  the  spacial  and  frequency  domains  starts  with  the  principle  of  duality  between  the  two — 
multiplication  in  the  time  domain  is  convolution  in  the  frequency  domain,  and  vice  versa.  Other 
relationships  between  the  two  domains  will  briefly  be  discussed  (Parseval's  Theorem,  etc.). 

3.1  Correlation 

Correlation  is  the  first  step  in  understanding  convolution.  An  excellent  reference  is  [15]  from  which  the 
following  definition  is  obtained  and  the  first  four  pages  are  included  in  Appendix  A. 3. 

i=-N 


Here  F  is  considered  of  odd  length.  Consider  it  a  sliding  window  operation — overlapping  the  sequences 
and  then  multiplying  and  adding.  Additional  coefficients  may  be  involved — to  perform  functions  such  as 
averaging,  etc.  When  F  does  not  completely  overlap  I  operations  are  undefined  for  values  outside  of  the 
boundaries  of  I.  Several  options  are  available. 

•  pad  with  zeros  (MATLAB  default) 

•  pad  with  1  st/last  values 

•  cyclically  repeated 

Computation  can  continue  when  undefined  values  are  resolved  with  one  of  these  methods. 

3.2  Convolution 

Convolution  is  the  dual  operation  to  multiplication  in  the  spacial  and  frequency  domains.  It  is  the  same  as 
correlation,  except  that  the  filter  is  flipped  first  [15]  (e.g.  F=(2,8,9)  ->  (9,8,2))  and  defined  as 

F*/(.r)=  f>(/)/(.r-0 

i=-N 

Another  useful  reference  on  convolution  is  [16]. 

3.3  Parseval’s  Theorem 

One  of  the  other  relationships  between  the  spacial  and  frequency  domain  is  Parseval’s  Theorem  [17],  It 
is  defined  as 


S  r*wr = 

which  shows  the  equivalence  between  energy  in  the  signal  regardless  of  the  domain.  It  is  used  to  check 
computational  results  to  ensure  the  principle  is  preserved. 

3.4  Discrete  Fourier  Transform 

The  Discrete  Fourier  Transform  (DFT)  is  a  mathematical  tool  [18]  used  to  change  data  representations 
between  the  spacial  and  frequency  domains.  A  variant  that  allows  faster  compute  times  is  called  the  Fast 
Fourier  Transform  (FFT)  but  performs  the  equivalent  function.  The  inverse  DFT  of  convolution  results  in 
the  frequency  domain  should  match  multiplication  results  in  the  spacial  domain.  This  method  will  be  used 
to  verify  frequency  domain  calculations. 


1  r¥— L 

—  T 

y 


3.4.1  Scaling 

From  [18]  the  following  is  informative  regarding  the  scale  factor  applied  to  the  DFT: 
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The  normalization  factor  multiplying  the  DFT  and  IDFT  (here  1  and  1/N)  and  the  signs  of 
the  exponents  are  merely  conventions,  and  differ  in  some  treatments.  The  only 
requirements  of  these  conventions  are  that  the  DFT  and  IDFT  have  opposite-sign 
exponents  and  that  the  product  of  their  normalization  factors  be  1/N.  A  normalization  of 


for  both  the  DFT  and  IDFT,  for  instance,  makes  the  transforms  unitary. 


Note  that  MATLAB  uses  the  opposite  scaling  mentioned  in  the  first  sentence.  From  [19]  elaboration  on 
scaling  is  given: 


Question  3.7.  FFTW  gives  results  different  from  my  old  FFT. 


People  follow  many  different  conventions  for  the  DFT,  and  you  should  be  sure  to  know 
the  ones  that  we  use  (described  in  the  FFTW  manual).  In  particular,  you  should  be  aware 
that  the  FFTW_FORWARD/FFTW_BACKWARD  directions  correspond  to  signs  of  -1/+1 
in  the  exponent  of  the  DFT  definition.  (Numerical  Recipes  uses  the  opposite  convention.). 
You  should  also  know  that  we  compute  an  unnormalized  transform.  In  contrast,  Matlab  is 
an  example  of  program  that  computes  a  normalized  transform.  See  Q3.10  'Why  does 
your  inverse  transform  return  a  scaled  result?'. 


Finally,  note  that  floating-point  arithmetic  is  not  exact,  so  different  FFT  algorithms  will  give 
slightly  different  results  (on  the  order  of  the  numerical  accuracy;  typically  a  fractional 
difference  of  1e-15  or  so  in  double  precision). 


Question  3.10.  Why  does  your  inverse  transform  return  a  scaled  result? 

Computing  the  forward  transform  followed  by  the  backward  transform  (or  vice  versa) 
yields  the  original  array  scaled  by  the  size  of  the  array.  (For  multi-dimensional  transforms, 
the  size  of  the  array  is  the  product  of  the  dimensions.)  We  could,  instead,  have  chosen  a 
normalization  that  would  have  returned  the  unsealed  array.  Or,  to  accommodate  the 
many  conventions  in  this  matter,  the  transform  routines  could  have  accepted  a  "scale 
factor"  parameter.  We  did  not  do  this,  however,  for  two  reasons.  First,  we  didn't  want  to 
sacrifice  performance  in  the  common  case  where  the  scale  factor  is  1.  Second,  in  real 
applications  the  FFT  is  followed  or  preceded  by  some  computation  on  the  data,  into 
which  the  scale  factor  can  typically  be  absorbed  at  little  or  no  cost. 

3.5  Fast  Fourier  Transform  (FFT)  Proper  Scaling  Test 

Figure  4  shows  the  proper  scaling  of  the  FFT  in  MATLAB  to  match  the  input  signal.  From  the  three  term 
sinusoidal  combination  the  frequency  domain  peak  values  match  those  defined  in  the  time  domain  signal. 


%  MATLAB  example 

Fsl  =  1000;  %  F=l/T  =l/le-3sec  =  1/msec 

L  =  1000; 

t=  (0:1/ Fs 1:5000-1) ;  %  5  sec  @  Fs=1000 

y  =  0.3*sin(2*pi*35*t)+0.7  *sin  (2*pi*50*t)  +  sin (2*pi*120*t) ; 
tNFFT  =  2 Anextpow2 (L) ; 
tY  =  f ft (y, tNFFT) /tNFFT; 

tf  =  Fsl/2*!3Lnspace  (0, 1,  tNFFT/2+1)  ;  %  one  sided  frequency 

subplot (2,1,1) ,plot(1000*t (1:250) ,y (1:250) ) 
subplot (2,1,2) , plot (tf, 2*abs (tY ( 1 : tNFFT/2+1 ) ) ) 
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FFT  Example  -  0.3cos(35Hz)+0.7cos(50Hz)+1  .Ocos(1 20Hz) 


Figure  4  -  FFT  Example  with  Proper  Scaling 


4.0  FREQUENCY  DOMAIN 

Application  of  Section  2.0  using  Section  3.0  will  be  demonstrated  in  Section  4.0. 

4.1  Root-Mean-Square  (RMS)  Terrain  Roughness  Metric 

The  calculation  of  the  RMS  terrain  roughness  metric  in  the  frequency  domain  is  based  on  the  duality  of 
convolution  (multiplication  in  the  spacial  domain).  The  approach  here  is  to  provide  equivalent  frequency 
domain  operations  and  examples  that  show  the  same  results  as  in  the  spacial  domain  with  a  plot  of 
results  and  description  of  the  MATLAB  script  used  to  generate  it.  It  should  be  noted  that  index  values  in 
MATLAB  are  one  (1)  based  and  not  zero  (0)  based. 

4.2  Munson  Gravel 

The  following  examples  apply  to  the  profile  known  as  Munson  Gravel. 

4.2.1  Spacial  Domain 

As  described  in  Section  2.3  the  time  domain  is  the  starting  point.  Figure  5  is  based  on  time  domain 
operations  (see  Appendix  A.1).  Notice  that  the  mean  and  detrended  profile  lose  one  half  the  length  of  the 
filter  at  the  beginning  and  end. 


y=dlmread ( ' MUNGL . FIL ' ) ;  %  read  in  the  data 

x  =  y(:,2)*12;  %  convert  feet  to  inches  or  consistent  units 

tmp  =  size (x)  ; 

lenx  =  tmp(l);  %  get  the  length  (#  of  points) 

%  e- (N*A/lambda) 

N  =  [0:1:120]  1 ;  %  1-side  =  30  ft  =  360  in  @  3  in  samples  =  120  (360/3) 

A  =  3;  %  3  in  samples 

lambda  =  10*12;  %  10  ft  *  12  =  120  in 

%  -N*A/lambda  =  -120*3/120  =  -3  @  last  sample 
[ fdet , rms , fmean, xmean] =rmswes (x,  N,  A, lambda) ; 

ttl  =  (0 : lenx-1)  ' *  (1/3) ;  %  displacement,  3  in  samples 

plot  (ttl , x,  ' r ' ) ; 
hold  on 

plot (ttl (121 : end-120) , fmean, '  g*  )  ;  %  account  for  filter  beg  and  end 

plot (ttl (121 : end-120) , fdet, 'b ') ;  %  account  for  filter  beg  and  end 


Munson  Gravel 


Figure  5  -  Spacial  Domain  RMS  Terrain  Roughness  Metric  Results 

Figure  6  shows  a  zoomed  view  of  Figure  5  for  the  terrain  profile  and  trend  based  on  the  filter  defined  in 
Section  2.2.  Notice  the  missing  samples  in  the  trend  if  the  original  data  is  not  extended  by  the  length  of 
the  filter  to  avoid  loss  of  data. 


plot (ttl , x, ' r ' ) ; 
hold  on 

plot (ttl ( 121 : end-120 ) , fmean, ' g 1 ) ; 
axis ( [1 . 444e4, 1 . 455e4 , 1732 , 1737 ] ) 
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Munson  Gravel 


Figure  6  -  Zoomed  Terrain  Profile  and  Trend 
Figure  7  shows  a  zoomed  view  of  the  detrended  terrain  profile. 


plot  (ttl,x,  ' r ' ) ; 
hold  on 

plot (ttl (121 : end- 120) ,  fmean, '  g' ) ; 
axis ( [1 . 444e4, 1 . 455e4 , 1732 , 1737 ] ) 
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Munson  Gravel 


Figure  7  -  Detrended  Terrain  Profile 


4.2.2  Spacial  Domain  Filters 

Figure  8  shows  normalized  2-sided  exponential  and  rectangular  filters.  Notice  the  overlap  at  fexp(lenN). 


%  2-sided  exponential 
tmp  =  size (N) ; 
lenN  =  tmp(l); 

fexp  =  zeros (lenN*2-l , 1 ) ;  %  overlap  at  zero 

fexp(l:lenN)  =  exp (-N (length (N) : -1 : 1 ) *A/lambda) ; 
fexp (lenN : end)  =  exp (-N*A/lambda) ; 

fexp (lenN)  =  2;  %  double  count  zero  position 

fexp  =  fexp/sum ( fexp) ; 
tmp  =  size (fexp); 
lenfexp  =  tmp(l); 

%  rect 

rect  =  ones (lenN*2-l , 1 ) ; 
rect  =  rect/sum (rect) ; 
tmp  =  size  (rect); 
lenrect  =  tmp(l); 

subplot  (2, 1 , 1 ) , plot (fexp,  ' r ' ) 
subplot  (2,1,2) , plot (rect ,  ' g ' ) 
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Spacial  Windows  -  Normalized 
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Figure  8  -  Spacial  Domain  Filters 


4.2.3  Frequency  Domain  Filters 

Figure  9  shows  the  frequency  domain  version  of  the  spacial  domain  filters.  In  real  frequency  terms, 
accounting  for  negative  frequencies,  only  one  half  of  the  length  (NFFT)  is  used  with  double  the  absolute 
value  of  the  magnitude. 


%  filters  in  frequency  domain  -  FFT 

NFFT  =  2Anextpow2 (lenx) ;  %  ! ! !use  all  with  this  so  can  mul  in  freq  domain 

Ffexp  =  f ft ( fexp, NFFT) /NFFT;  %  !!! scale  by  len  to  get  right  mag  (see  sin  test) 
Frect  =  fft (rect,NFFT) /NFFT; 

%  Fs  =  1/T 

Fs=(l/3);  %  4  samples  in  12  in  =  4/12  =  1/3 

%  linspace (xl , x2 , N)  =  N  pts  between  xl  and  x2 
f  =  Fs/2*linspace (0, 1, NFFT/2+1) ; 

subplot (2,1,1) , plot (f , 2*abs (Ffexp (1 : NFFT/ 2+1) )  , ' r' ) ;  %  !!!  2*abs  /NFFT/2+1 

subplot (2,1,2) , plot (f , 2*abs (Frect (1 : NFFT/2+1) ) , ' g' ) ; 
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Frequency  Windows  (FFT) 
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Figure  9  -  Frequency  Domain  Filters 

4.2.3. 1  Spacial  Domain  vs.  Frequency  Domain  Filter  Differences 

Based  on  proper  scaling  demonstrated  in  Section  3.5,  differences  in  the  spacial  domain  and  inverse  FFT 
of  the  frequency  domain  defined  filters  are  shown  in  Figure  10.  Note  they  are  the  same  except  for  round 
off  error. 


%  filter  differences  -  FFT/inverse  FFT 

ffexp  =  NFFT*if f t (Ffexp, NFFT) ;  %!!!  scale  by  length 

frect  =  NFFT*ifft (Frect, NFFT) ; 

subplot (2,1,1) , plot (fexp-ffexp (1: lenfexp) , ' r ’ ) ;  %  lenfexp  =  lenN*2-l 

subplot (2,1,2) , plot (rect- frect ( 1 : lenrect) , ' b ' ) ; 
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Spacial/Frequency(FFT/iFFT)  Window  Differences 
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Figure  10  -  Spacial  versus  Frequency  Domain  Filter  Differences 

4.2.4  Frequency  Domain  Operations  -  FFT/iFFT 

Figure  11  shows  the  error  between  the  original  terrain  profile  and  the  FFT/iFFT  of  it — indicating  the 
operation  gives  back  the  original  with  no  error. 


%  FFT  x 

Fx  =  f ft (x, NFFT) /NFFT; 

tt6= ( 0 : lenx-1 ) * (1/3) ; 

ft6  =  Fs/2*linspace (0, l,NFFT/2+l) ; 

%  inverse  FFT  x 
Fcx  =  Fx; 

f Fcx  =  NFFT*ifft (Fcx, NFFT) ; 

plot  (tt6, x,  ' r ' ) ; 
hold  on 

plot (tt6, f Fcx ( 1 : lenx) , ' g* ') ;  %  !!!  NFFT  ->  true  length 
plot (tt6, x-fFcx (1 : lenx) , '  b ' ) ; 
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Munson  Gravel 


Figure  11  -  Frequency  Domain  Operations  -  FFT/iFFT  Comparison 

Figure  12  zooms  in  on  the  end  of  the  original  and  FFT/iFFT  terrain  profile  to  show  both  are  equivalent, 
validating  the  representation  of  the  terrain  profile  in  the  spacial  or  frequency  domain. 
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Munson  Gravel 
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Figure  12  -  Zoomed  FFT/iFFT  Comparison 

4.2.5  Spacial  vs.  Frequency  Domain  Mean  Comparison 

Figure  13  shows  the  comparison  of  the  terrain  profile  trend  (or  mean)  where  it  was  filtered  in  the  spacial 
and  frequency  domain  showing  the  equivalence  of  multiplication  in  the  spacial  domain  and  convolution  in 
the  frequency  domain. 


%  exp  filter 
H  =  Ffexp; 
tmp  =  size  (H) ; 
lenH  =  tmp(l); 

%  inverse  x,  mean,  mydet,  detrended 
Fcx  =  Fx.*H; 

f Fcx  =  NFFT*NFFT*ifft (Fcx,NFFT) ;  %  !!!  scale  by  time  length 

plot ( fmean, ' r ' ) 
hold  on 

%  ! ! !  proper  indexes  for  y  cutoffs  (y  only  when  exp  fully  in  orig  y) 
plot (abs (fFcx (lenfexp : lenx+ (lenfexp-1) - (lenfexp-1) ) ) , ' g* ' ) ; 
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Munson  Gravel  -  Mean 


A  zoomed  version  is  shown  in  Figure  14  for  the  end  of  the  terrain  profile. 


Munson  Gravel  -  Mean 


Figure  14  -  Zoomed  Mean  Comparison 

4.2.6  Spacial  vs.  Frequency  Domain  Detrended  Comparison 

Figure  15  shows  the  detrended  terrain  profile  based  on  subtracting  the  trend  obtained  from  the  frequency 
domain. 


%  detrended 

%  assumes  lenexp  is  odd,  take  out  mean  from  FFT/ . *H/iFFT 

mydet  =  x(  (lenfexp+1 ) /2  :  lenx  -  (lenfexp+1) /2  +  1  )  -  abs (  f Fcx (  lenfexp  :  lenx  +  (lenfexp- 

1)  -  (lenfexp-1)  )  ) ; 

mm=me  an ( myde  t ) 
mydet  =  mydet  -  mm; 

plot (fdet, ' r ' ) ; 
hold  on 

plot (  mydet,' g* 1 ) 
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Munson  Gravel  -  Detrended 


0  0.5  1  1.5  2  2.5  3  3.5  4  4.5 

Distance  [in]  4 

x  10 

Figure  15  -  Spacial  versus  Frequency  Domain  Detrended  Comparison 


A  zoomed  version  of  the  end  is  shown  in  Figure  16.  Note  the  agreement  in  results. 


Munson  Gravel  -  Detrended 
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Figure  16  -  Zoomed  Detrended  Comparison 

4.2.7  Spacial  vs.  Frequency  Domain  Mean  and  Detrended  Differences 

Figure  17  shows  the  difference  between  the  spacial  and  frequency  domain  for  the  mean  and  detrended 

profile.  Note  they  are  the  same  except  for  round  off  error — about  machine  precision. 


%  mean  error 

%  ! ! !  proper  indexes  for  y  cutoffs  (y  only  when  exp  fully  in  orig  y) 

subplot ( 2 , 1, 1) ,plot (fmean  -  abs (fFcx (lenfexp: lenx+ (lenfexp-1) - (lenfexp-1) ) ) ,  ' r ' ) ; 

%  detrended  error 

subplot (2,1,2) , plot ( fdet-mydet , ' r ' ) ; 
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Figure  17  -  Spacial  versus  Frequency  Domain  Mean/Detrended  Differences 

4.2.8  Equivalent  Measures 

Results  output  from  the  script  below  show  numerical  comparisons  between  the  spacial  domain  RMS 
function,  standard  deviation  of  the  detrended  signal,  square  root  of  the  variance  of  the  detrended  signal, 
and  each  side  of  Parseval’s  Theorem  outlined  in  Section  3.3  to  be  equivalent. 


%  PSD 

%  std(x,l)  =  sqrt ( sum (x . A2 ) /length (x) ) 

tmp  =  size (mydet) ; 
lenmydet  =  tmp(l) 

Fmydet  =  f ft (mydet , NFFT) /NFFT; 

Pxx  =  Fmydet. *conj (Fmydet) ; 

! rm  tlO-d 
diary  tlO-d 
'  spacial/freq  results' 

'  rms  ' 

' std (mydet, 1 ) ' 

' sqrt (var (mydet, 1 ) ) ' 

' sqrt (  sum (mydet . A2 ) / length (mydet) ) ' 

'sqrt(  sum (  Pxx  ) ^length (mydet) /NFFT) ' 

[rms, std (mydet, 1) , sqrt (var (mydet, 1) ), sqrt (  sum (mydet . A2 ) /length (mydet )), sqrt (  sum  (  Pxx 

) *NFFT/ length (mydet) ) ] 

diary 
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ans  = 

spacial/freq  results 

ans  = 

rms 


ans 

s  t  d (mydet , 1 ) 

ans  = 

sqrt (var (mydet, 1) ) 

ans  = 

sqrt (  sum (mydet .  A2 )  / length (mydet) ) 

ans  = 

sqrt (  sum (  Pxx  ) *length (mydet ) /NFFT) 

ans  = 

0.3364  0.3364  0.3364  0.3364  0.3364 


Computation  of  the  Root-Mean-Square  (RMS)  Terrain  Roughness  Metric  (Section  2.1)  is  valid  in  either 
the  spacial  or  frequency  domain. 

4.3  Perryman  3 

The  following  examples  apply  to  the  profile  known  as  Perryman  3.  Differences  from  Section  4.2  are  only 
in  the  file  read  in  and  label  names. 

4.3.1  Spacial  Domain 

As  described  in  Section  2.3  the  time  domain  is  the  starting  point.  Figure  18  is  based  on  time  domain 
operations. 
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Perryman  3 
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Figure  18  -  Spacial  Domain  RMS  Terrain  Roughness  Metric  Results 

Figure  19  shows  a  zoomed  view  of  Figure  18  for  the  terrain  profile  and  trend  based  on  the  filter  defined  in 
Section  2.2.  Notice  the  missing  samples  in  the  trend  if  the  original  data  is  not  extended  by  the  length  of 
the  filter  to  avoid  loss  of  data. 
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Perryman  3 
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Figure  19  -  Zoomed  Terrain  Profile  and  Trend 
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Figure  20  shows  a  zoomed  view  of  the  detrended  terrain  profile. 


Perryman  3 


Figure  20  -  Detrended  Terrain  Profile 


4.3.2  Frequency  Domain  Operations  -  FFT/iFFT 

Figure  21  shows  the  error  between  the  original  terrain  profile  and  the  FFT/iFFT  of  it — indicating  the 
operation  gives  back  the  original  with  no  error. 
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Perryman  3 


Figure  21  -  Frequency  Domain  Operations  -  FFT/iFFT  Comparison 

Figure  22  zooms  in  on  the  end  of  the  original  and  FFT/iFFT  terrain  profile  to  show  both  are  equivalent, 
validating  the  representation  of  the  terrain  profile  in  the  spacial  or  frequency  domain. 
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Perryman  3 


Figure  22  -  Zoomed  FFT/iFFT  Comparison 

4.3.3  Spacial  vs.  Frequency  Domain  Mean  Comparison 

Figure  23  shows  the  comparison  of  the  terrain  profile  trend  (or  mean)  where  it  was  filtered  in  the  spacial 
and  frequency  domain  showing  the  equivalence  of  multiplication  in  the  spacial  domain  and  convolution  in 
the  frequency  domain. 
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Perryman  3  -  Mean 
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Figure  23  -  Spacial  versus  Frequency  Domain  Mean  Comparison 
A  zoomed  version  is  shown  in  Figure  24  for  the  end  of  the  terrain  profile. 
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Perryman  3  -  Mean 


Figure  24  -  Zoomed  Mean  Comparison 

4.3.4  Spacial  vs.  Frequency  Domain  Detrend  Comparison 

Figure  25  shows  the  detrended  terrain  profile  based  on  subtracting  the  trend  from  the  frequency  domain. 
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Figure  25  -  Spacial  versus  Frequency  Domain  Detrended  Comparison 


A  zoomed  version  of  the  end  is  shown  in  Figure  26.  Note  the  agreement  in  results. 


Perryman  3  -  Detrended 


Figure  26  -  Zoomed  Detrended  Comparison 

4.3.5  Spacial  vs.  Frequency  Domain  Mean  and  Detrended  Differences 

Figure  27  shows  the  difference  between  the  spacial  and  frequency  domain  for  the  mean  and  detrended 

profile.  Note  they  are  the  same  except  for  round  off  error. 
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Figure  27  -  Spacial  versus  Frequency  Domain  Mean/Detrended  Differences 

4.3.6  Equivalent  Measures 

Results  output  from  the  script  in  Section  4.2.8  for  this  terrain  profile  show  numerical  comparisons  between 
the  spacial  domain  RMS  function,  standard  deviation  of  the  detrended  signal,  square  root  of  the  variance 
of  the  detrended  signal,  and  each  side  of  Parseval’s  Theorem  outlined  in  Section  3.3  to  be  equivalent. 
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ans  = 

spacial/freq  results 

ans  = 

rms 


ans 

s  t  d (mydet , 1 ) 

ans  = 

sqrt (var (mydet, 1) ) 

ans  = 

sqrt (  sum (mydet . A2 ) / length (mydet) ) 

ans  = 

sqrt (  sum (  Pxx  ) *length (mydet ) /NFFT) 

ans  = 

3.1228  3.1228  3.1228  3.1228  3.1228 


Again,  it  is  demonstrated  (with  a  completely  different  terrain  profile)  that  the  computation  of  the  Root- 
Mean-Square  (RMS)  Terrain  Roughness  Metric  (Section  2.1)  is  valid  in  either  the  spacial  or  frequency 
domain. 

4.4  Test  Case  Calculations 

An  example  test  case  is  included  in  Appendix  A. 5  where  corresponding  script  functions  and  output  can  be 
compared  for  a  simple  signal 

x  =  [1  3  -10  5  8  9  22  2  3  10  10  12  11  14  13  9  17] 


and  filter 


h  =  [123832  1]’; 

to  demonstrate  how  MATLAB  functions  (see  Appendix  A.2)  perform. 
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APPENDIX  A.1  -  RMS  (Exponential  Weighted  Filter) 


%function  [fdet, rms, fmean, xmean] =rms (x,N,  A, lambda) 

%  calculate  RMS  exp. weighted  average  value  of  terrain 

%inputs : 

%  x  -  array  to  filter 

%  N  -  number  of  samples  vector  eg.  [012345] 

%  A  -  sample  spacing 
%  lambda  -  exp  weighting  constant 

^returns : 

%  fdet  -  detrended  terrain 
%  rms  -  rms  value 
%  fmean  -  trend 
%  xmean  -  mean  of  fdet 

function  [fdet, rms , fmean, xmean] =rms (x, N, A, lambda) ; 

%  can  extend  input  by  length  of  filter  so  do  not  lose  data 

%x  =  [x ( 1 ) *zeros ( 1 , length (N) -1 )  ,  x',  x (end) *zeros ( 1 , length (N) -1 )]' ; 

n  =  length (x) ; 

s  =  length (N) ; 

e  =  (n+1)  -  s; 

detl  =  exp (-N*A/ lambda ) ; 

d  =  sum (detl ) *2 ; 

fmean  =  x; 

fdet  =  x; 

for  i  =  s:e, 

fmean(i)  =  sum (x (i : -1 : i-s+1 ) ,*detl  +  x (i : i+s-1) . *detl) /d; 
end 

fdet(s:e)  =  x(s:e)  -  fmean (s : e); 
xmean  =  mean (fdet (s : e) ) ; 
fdet(s:e)  =  fdet(s:e)-  xmean; 

%  for  best  unbiased  estimator  use  std(x,0)  to  divide  by  N-l  instead  of  N 
rms  =  std (fdet (s : e) , 1) ; 

fmean  =  fmean (s:e); 
fdet  =  fdet  (s : e) ; 
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APPENDIX  A.2  -  MATLAB  Function  Definitions 


dlmread  (http://www.mathworks.com/help/matlab/ref/dlmread.html) 
size  (http://www.mathworks.com/help/matlab/ref/size.html) 
nextpow2  (http://www.mathworks.com/help/matlab/ref/nextpow2.html) 
conv  (http://www.mathworks.com/help/matlab/ref/conv.html) 
fft  (http://www.mathworks.com/help/matlab/ref/fft.html) 

ifft  (http://www.mathworks.com/help/matlab/ref/ifft.html) 

The  functions  Y  =  fft(x)  and  y  =  ifft(X)  implement  the  transform  and  inverse 
transform  pair  given  for  vectors  of  length  N  by: 

Xik)  =  £  .v(/)£^-1:ii:*-1:i 
j=  1 

jV 

xtj)  =  (1  /  N ) £  * **-1 ' 

*=1 


where 


0* 


JV 


is  an  A/th  root  of  unity. 


var  (http://www.mathworks.com/help/matlab/ref/var.html) 

std  (http://www.mathworks.com/help/matlab/ref/std.html) 
conj  (http://www.mathworks.com/help/matlab/ref/coni.html) 
exp  (http://www.mathworks.com/help/matlab/ref/exp.html) 
linspace  (http://www.mathworks.com/help/matlab/ref/linspace.html) 
psd  (http://www.mathworks.com/help/signal/ref/spectrum.html) 
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APPENDIX  A.3  -  Reference  [15]  -  First  4  pages 


Correlation  and  Convolution 

Class  Notes  for  CMSC  426.  Fall  2005 
David  Jacobs 

Introduction 

Correlation  and  Convolution  are  basic  operations  that  we  will  perform  to  extract 
information  from  images.  They  are  in  some  sense  the  simplest  operations  that  we  can 
perform  on  an  image,  but  they  are  extremely  useful.  Moreover,  because  they  are  simple, 
they  can  be  analyzed  and  understood  very  well  and  they  are  also  easy  to  implement  and 
can  be  computed  very  efficiently.  Our  main  goal  is  to  understand  exactly  what 
correlation  and  convolution  do.  and  why  they  are  useful.  We  will  also  touch  on  some  of 
their  interesting  theoretical  properties;  though  developing  a  full  understanding  of  them 
would  take  more  time  than  we  have. 

These  operations  have  two  key  features:  they  are  shift-invariant.  and  they  are  linear. 
Shift-invariant  means  that  we  perform  the  same  operation  at  every  point  in  the  image. 
Linear  means  that  this  operation  is  linear,  that  is,  we  replace  every  pixel  with  a  linear 
combination  of  its  neighbors.  These  two  properties  make  these  operations  very  simple; 
it’s  simpler  if  we  do  the  same  thing  everywhere,  and  linear  operations  are  always  the 
simplest  ones. 

We  will  first  consider  the  easiest  versions  of  these  operations,  and  then  generalize.  Well 
make  things  easier  in  a  couple  of  ways.  First,  convolution  and  correlation  are  almost 
identical  operations,  but  students  seem  to  find  convolution  more  confusing.  So  we  will 
begin  by  only  speaking  of  correlation,  and  then  Later  describe  convolution.  Second,  we 
wiU  start  out  by  discussing  ID  images.  We  can  think  of  a  ID  image  as  just  a  single  row 
of  pixels.  Sometimes  things  become  much  more  complicated  in  2D  than  ID.  but  luckily, 
correlation  and  convolution  do  not  change  much  with  the  dimension  of  the  image,  so 
understanding  things  in  ID  will  help  a  lot.  Also,  later  we  will  find  that  in  some  cases  it  is 
enlightening  to  think  of  an  image  as  a  continuous  function,  but  we  will  begin  by 
considering  an  image  as  discrete,  meaning  as  composed  of  a  collection  of  pixels. 

Notation 

We  will  use  uppercase  letters  such  as  /  and  J to  denote  an  image.  An  image  may  be 
either  2D  (as  it  is  in  real  life)  or  ID.  We  wiU  use  lowercase  letters,  like  i  andj  to  denote 
indices,  or  positions,  in  the  image.  When  we  index  into  an  image,  we  will  use  the  same 
conventions  as  MatJab.  First,  that  means  that  the  first  element  of  an  image  is  indicated  by 
1  (not  0,  as  in  Java,  say).  So  if  J  is  a  ID  image,  1(1)  is  its  first  element.  Second,  for  2D 
images  we  give  first  the  row.  then  the  column.  So  1(3 f  6)  is  the  pixel  in  the  third  row  of 
the  image,  and  the  sixth  column. 


An  Example 
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One  of  the  simplest  operations  that  we  can  perform  with  correlation  is  local  averaging. 
As  we  will  see,  this  is  also  an  extremely  useful  operation.  Let's  consider  a  simple 
averaging  operation,  in  which  we  replace  every  pixel  in  a  ID  image  by  the  average  of 
that  pixel  and  its  two  neighbors.  Suppose  we  have  an  image  /  equal  to: 


5 

4 

2 

3 

7 

4 

6 

5 

3 

6 

Averaging  is  an  operation  that  takes  .an  image  as  input,  and  produces  a  new  image  as 
output.  When  we  average  rhe  fourth  pixel,  for  example,  we  replace  the  value  3  with  the 
average  of  2a  3,  and  7.  That  is,  if  we  call  the  new  image  that  we  produce  J  we  can  write: 
J(4)  =  (T(3)+I(4)+I(5))/3  =  (2+3+7)/3  =  4.  Or,  for  example,  we  also  get:  J(3)  = 
(I(2)-l(3)-I(4))/3  =  (4-2-3)73  =  3.  Notice  that  every  pixel  in  the  new  image  depends 
on  the  pixels  in  the  old  image.  A  possible  error  is  to  use  J(3)  when  calculating  J(4). 

Don't  do  this;  1(4)  should  only  depend  on  1(3),  1(4)  and  1(5).  Averaging  like  this  is  shift- 
invariant.  because  wre  perform  the  same  operation  at  every  pixel.  Every  new  pixel  is  the 
average  of  itself  and  its  two  neighbors.  Averaging  is  linear  because  every  new  pixel  is  a 
linear  combination  of  the  old  pixels.  This  means  that  we  scale  the  old  pixels  (in  this 
case,  we  multiply  all  the  neighboring  pixels  by  1/3)  and  add  them  up.  This  example 
illustrates  another  property  of  all  correlation  and  convolution  that  we  will  consider.  The 
output  image  at  a  pixel  is  based  on  only  a  small  neighborhood  of  pixels  around  it  in  the 
input  image.  In  this  case  the  neighborhood  contains  only  three  pixels.  Sometimes  we 
will  use  slightly  larger  neighborhoods,  but  generally  they  will  not  be  too  big. 

Boundaries:  We  still  haven't  fully  described  correlation,  because  we  haven't  said  what 
to  do  at  the  boundaries  of  the  image.  What  is  J(l)l  There  is  no  pixel  on  its  left  to  include 
in  the  average,  ie.r  1(0)  is  not  defined.  There  are  four  common  ways  of  dealing  with  this 
issue. 
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In  the  first  method  of  handling  boundaries,  the  original  image  is  padded  with  zeros 
(in  red  italics). 


The  first  way  is  to  imagine  that  /  is  part  of  an  infinitely  long  image  which  is  zero 
everywhere  except  where  we  have  specified.  In  that  case,  we  have  1(0}  =  ft,  and  we  can 
say:  J(l)  =  (1(0)  - 1(1)  + 1(2})/3  =  (0  +  5  +  4)/3  =  3.  Similarly,  we  have:  JflQ)  = 

(1(9)  +1(1Q)+I(1 l))/3  =(3  +  6+  0}/3  =  3. 
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In  the  second  method  of  handling  boundaries,  the  original  image  is  padded  with  the 
first  and  last  values  (in  red  italics). 


The  second  wray  is  to  also  imagine  that  1  is  part  of  an  infinite  image,  but  to  extend  it  using 
the  first  and  last  pixels  in  the  image.  In  our  example,  any  pixel  to  rhe  left  of  the  first  pixel 
in  /  would  have  the  value  5,  and  any  pixel  to  the  rielit  of  the  last  pixel  would  have  the 
value  6.  So  we  would  say:  J(I)  =  (1(0)  - 1(1)  - 1(2))73  =  (5  +  5  +  4)/3  =  4  2/3,  and 
1(10)  =  (1(9)-I(10)+1(U))73  =  (3-6-  6)73  =  5. 
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Ill  the  third  method  of  handling  boundaries,  the  original  image  is  repeated  cyclically 
(in  red  italics). 


Third,  we  can  imagine  the  image  as  being  like  a  circle,  so  that  the  pixel  values  repeat 
over  and  over  again.  The  pixel  to  the  left  of  the  first  pixel,  then,  would  be  the  last  pixel 
in  the  image.  That  is.  in  our  example,  we  would  define  1(0)  to  be  1(10).  Then  we  would 
have  J(l)  =  (1(0)  + 1(1)  + I(2))/3=  (1(10)  + 1(1 )- I(2))/3  =(6+5+  4)/3  =  5,  and  J(1 0) 
=  (lf9)-l(10)+l(ll))/3  =  (1(9)  -1(1 0)+l(l  ))/3  =  (3-6-  5)/3  =  4  2/3. 

Finally,  we  can  simply  say  that  the  image  is  undefined  beyond  the  values  that  we  have 
been  given.  In  that  case,  we  cannot  compute  any  average  that  uses  these  undefined 
values,  so  Jfl )  and  Jfl  0)  will  be  undefined,  and  J  will  be  smaller  than  /. 

These  four  methods  have  different  advantages  and  disadvantages  If  we  imagine  that  the 
image  we  are  using  is  j  ust  a  small  window  on  the  world,  and  we  want  to  use  values 
outside  the  boundary  that  are  most  similar  to  the  values  that  we  would  have  obtained  if 
we'd  taken  a  bigger  picture,  than  the  second  approach  probably  makes  the  most  sense. 
That  is,  if  we  had  to  guess  at  the  value  of  If  0),  even  though  we  can't  see  It,  the  value  we 
can  see  in  1(1)  is  probably  a  pretty  good  guess.  In  this  class,  unless  we  explicitly  state 
otherwise,  you  should  use  the  second  method  for  handling  boundaries. 

Correlation  as  a  Sliding,  Window  ed  Operation 

We're  now  going  to  look  at  the  same  averaging  operation  in  a  slightly  different  way 
which  is  more  graphical,  and  perhaps  more  intuitive  to  generalize,  In  averaging,  for  a 
specific  pixel  we  multiply  it  and  its  neighbors  by  1/3  each,  and  then  add  up  the  three 
resulting  numbers.  The  numbers  we  multiply,  (ITS.  IT.  1/3)  form  a  .filter.  This  particular 
filter  is  called  a  box  filter.  We  can  think  of  it  as  a  1x3  structure  that  we  slide  along  the 
image.  At  ea  ch  position,  we  multiply  each  number  of  the  filter  by  the  image  number  that 
lies  underneath  it,  and  add  these  all  up.  The  result  is  a  new  number  corresponding  to  the 
pixel  that  is  underneath  the  center  of  the  filter.  The  figure  below  shows  us  producing  Jfl) 
in  this  way. 
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To  produce  the  next  number  in  the  filtered  image,  we  slide  the  filter  over  a  pixel,  and 
perform  die  same  operation. 
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We  continue  doing  tins  until  we  have  produced  every  pixel  in  /  With  this  view  of 
correlation,  we  can  define  a  new  averaging  procedure  by  just  defining  a  new'  filter.  For 
example,  suppose  instead  of  averaging  a  pixel  with  its  immediate  neighbors,  we  want  to 
average  each  pixel  with  immediate  neighbors  and  their  immediate  neighbors.  We  can 
define  a  filter  as  (1/5, 1/5, 1/5,  1/5, 1/5).  Then  we  perform  the  same  operation  as  above, 
but  usinE  a  filter  that  is  five  pixels  wide.  The  first  pixel  in  the  resulting  image  will  then 
be:  JfU"=  ilf-l i  '5  -  1<0,  5  - 1(1  >  5  -  1C>.  3  - 1/3)  3/  =  1-1 -1-4. 3  -3  5  =  4  1  5. 


A  Mathematical  Definition  for  Correlation 

If  s  helpful  to  write  this  all  down  more  formally.  Suppose  F  is  a  correlation  filter.  It  wrill 
be  convenient  notationally  to  suppose  that  F  has  an  odd  number  of  elements,  so  we  can 
suppose  that  as  it  shifts,  its  center  is  right  on  top  of  an  element  of/  So  we  say  that  F  has 
2jV+/  elements,  and  that  these  are  indexed  from  -N to  N,  so  that  the  center  element  of  F  is 
F(O).  Then  we  can  write: 


N 

F°I(x)  =  y  F(i)I(x  +  i) 

z"=— jV 


where  the  circle  denotes  correlation.  With  this  notation,  we  can  define  a  simple  box  filter 

as: 


0 


for  /  = -1.0,1 
for  /  ^  -1,0,1 


Constructing  an  Filter  from  a  C  ontinuous  Function 

It  is  pretty  intuitive  what  a  reasonable  averaging  filter  should  look  like.  Now  we  want  to 
start  to  consider  more  general  strategies  for  constructing  filters.  It  commonly  occurs  that 
we  have  in  mind  a  continuous  function  that  would  make  a  good  filter,  .and  we  want  to 
come  up  with  a  discrete  filter  that  approximates  this  continuous  function.  Some  reasons 
for  thinking  of  filters  first  as  continuous  functions  will  be  given  when  we  talk  about  the 


41 


APPENDIX  A.4  -  Plot  Demonstration  Script 


%ti 

close  all 
clear  all 
fig  =  1 

y=dlmread ( ' MUNGL . FIL '  )  ; 

%y  =  dlmread ( ' 3courL . FIL ' )  ; 
x  =  y(:,2)*12;  %  convert  feet  to  inches 
tmp  =  size (x) ; 
lenx  =  tmp ( 1 ) ; 

%  e- (N*A/lambda) 

N  =  [0:1:120]';  %  1-side  =  30  ft  =  360  in  @  3  in  samples  =  120  (360/3) 
A  =  3;  %  3  in  samples 

lambda  =  10*12;  %  10  ft  *  12  =  120  in 
%  -N*A/lambda  =  -120*3/120  =  -3  @  last  sample 
% [ fdet , rms , fmean, xmean] =rmswes (x, N, A, lambda) ; 

[ fdet , rms , fmean, xmean] =rmswes (x, N, 3 , 10*12 ) ;  %  3  inch  samples  to  30  ft 

figure (fig) 
cl  f 

fig  =  fig+1 

ttl  =  ( 0 : lenx-1 )  '* (1/3) ; 
plot (ttl , x,  ' r ' )  ; 
hold  on 

plot (ttl (121 : end- 120) ,  fmean, 'g' ) ; 
plot  (ttl (121 : end- 120) , fdet,  'b ' ) ; 
title (' Munson  Gravel'); 

%title ( ' Perryman  3 ' ) ; 
ylabel (' Elevation  [in]'); 
xlabel (' Distance  [in]'); 

legend ( ' original ' , ' mean ' ,  ' detrended ' ,  ' Location ' ,  ' Northwest ' ) ; 
axis ( [-50,15000,-200,1800] ) 

print ( ' -dpng ' , ' tl-mun ' ) ; 

%print ( ' -dpng ' , ' tl-p3 ' ) ; 


figure (fig) 
cl  f 

fig  =  fig+1 

ttl  =  ( 0 : lenx-1 )  '* (1/3) ; 
plot (ttl , x,  ' r ' )  ; 
hold  on 

plot  (ttl (121 : end- 12  0) , fmean,  'g' ) ; 
title (' Munson  Gravel'); 

%title ( ' Perryman  3 ' )  ; 
ylabel (' Elevation  [in]'); 
xlabel (' Distance  [in]'); 

legend ( ' original ' , ' mean ' , ' Location ' , ' Northwest ' ) ; 
axis ( [1.444e4,1.455e4, 1732, 1737] ) 
print ( ' -dpng ' , ' tl-zml ' ) ; 

figure (fig) 
cl  f 

fig  =  fig+1 

ttl  =  ( 0 : lenx-1 )  '* (1/3) ; 

plot (ttl (121 : end- 120) , fdet, 'b ' ) ; 

title (' Munson  Gravel'); 

%title ( ' Perryman  3 ' ) ; 
ylabel (' Elevation  [in]'); 
xlabel (' Distance  [in]'); 

legend ( ' detrended ' , ' Location ' , ' Northwest ' ) ; 
axis ( [1.42e4,1.4  7e4,-4,3]  ) 
print ( ' -dpng ' , ' tl-zm2 ' ) ; 

%t2  .  m 


%  create  and  plot  filters 

%  2-sided  exponential 
tmp  =  size (N) ; 
lenN  =  tmp(l); 
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fexp  =  zeros (lenN*2-l , 1 ) ;  %  overlap  at  zero 

fexp(lilenN)  =  exp (-N (length (N) : -1 : 1 ) *A/lambda) ; 

fexp (lenN : end)  =  exp (-N*A/lambda) ; 

fexp(lenN)  =  2;  %  double  count  zero  position 

fexp  =  fexp/sum ( fexp) ; 

tmp  =  size (fexp); 

lenfexp  =  tmp(l); 

%  rect 

rect  =  ones (lenN*2-l , 1 ) ; 
rect  =  rect/sum (rect ) ; 
tmp  =  size  (rect) ; 
lenrect  =  tmp(l); 

figure (fig) 
cl  f 

fig  =  fig  +  1 

subplot (2,1,1) , plot (fexp, ' r ' ) 
ylabel ( ' exp ' ) 

title  (' Spacial  Windows  -  Normalized') 

axis (  [-1,252,-0.01,0.03]  )  ; 

xlabel (' Distance  [in]') 

subplot (2,1,2) , plot (rect, ' g ' ) 

ylabel ( ' rect ' ) ; 

xlabel ( ' Distance  [in]  '  ) 

axis ( [-1,252,-0.01,0.01] ) ; 

print ( ' -dpng ' , ' t2-tw ' ) ; 

%t3  .m 


%  plot  filters  FFT 
lenx 

NFFT  =  2Anextpow2 (lenx) ;  %  ! ! !use  all  with  this  so  can  mul  in  freq  domain 

lenfexp 

lenrect 

Ffexp  =  f ft ( fexp, NFFT) /NFFT;  %  !!! scale  by  len  to  get  right  mag  (see  sin  test) 
Frect  =  f ft (rect , NFFT) /NFFT; 

%  Fs  =  1/T 

Fs=(l/3);  %  4  samples  in  12  in  =  4/12  =  1/3 

%  linspace (xl , x2 , N)  =  N  pts  between  xl  and  x2 
f  =  Fs/2* linspace (0, 1, NFFT/ 2+1) ;  % 

figure (fig) 
cl  f 

fig  =  fig  +  1 

subplot (2, 1, 1) , plot (f, 2*abs (Ffexp (1 :NFFT/2+l) ),' r ') ;  %  !!!  2*abs  /NFFT/2+1 

title (' Frequency  Windows  (FFT) '); 

ylabel ( ' exp ' ) 

xlabel (' Frequency  [Hz]') 

axis ([0,0.05,0, 4e-5] ) ; 

subplot (2,1,2) , plot (f , 2*abs (Frect (1 : NFFT/2+1) ) , ' g' ) ; 

ylabel ( ' rect ' ) 

xlabel (' Frequency  [Hz]') 

axis ([0,0.05,0, 4e-5]  )  ; 

print ( ' -dpng ' , ' t3-f f t ' ) ; 

%t4  .  m 


%  plot  cos  test  FFT 

%  MATLAB  example 
figure (fig) 
cl  f 

fig  =  fig  +  1 

Fsl  =  1000  %  F=l/T  =l/le-3sec  =  1/msec 
L  =  1000; 

t=  (0:1/ Fs 1:5000-1) ;  %  5  sec  @  Fs=1000 


y  =  0 . 3*sin (2*pi*35*t) +0 . 7  *sin (2*pi*50*t)  +  sin (2*pi*120*t) ; 
tNFFT  =  2 Anextpow2 (L) ; 
tY  =  f ft (y, tNFFT) /tNFFT; 

tf  =  Fsl/2*linspace (0, 1, tNFFT/2+1) ;  %  one  sided  frequency 

subplot (2,1,1) , plot (1000*t (1:250) ,y (1:250) ) 

title('FFT  Example  -  0 . 3cos (35Hz) +0 . 7cos (50Hz) +1 . Ocos (120Hz) ' ) ; 
ylabel ( ' Signal ' ) 
xlabel ( 'Time  [msec] ' ) ; 

subplot (2,1,2) , plot (tf , 2*abs (tY ( 1 : tNFFT/2+1 ) ) ) 
xlabel (' Frequency  [Hz]'); 
ylabel ( ' FFT  (2* | Y (f ) | ) ' ) ; 
print ( ' -dpng ' , ' t4-cos ' ) ; 

%t5  .  m 


%  plot  filter  inverse  FFT 

ffexp  =  NFFT*if f t (Ffexp, NFFT) ;  %  !!!  scale  by  length 
frect  =  NFFT*ifft (Frect,NFFT) ; 

figure (fig) 
cl  f 

fig  =  fig  +  1 

subplot (2 , 1 , 1 ), plot ( fexp-ffexp ( 1 : lenfexp) ,' r ') ;  %  lenfexp  =  lenN*2-l 
title (' Spacial/Frequency (FFT/iFFT)  Window  Differences'); 
ylabel ( ' exp ' ) ; 

subplot (2,1,2) , plot (rect- frect ( 1 : lenrect ) , ' b ' ) ; 

ylabel  (  ' rect ' ) 

print ( ' -dpng ' , ' t5-wd ' ) ; 

%t6 . m 


%  plot  FFT  x 

Fx  =  fft (x, NFFT) /NFFT; 

tt6= ( 0 : lenx-1 ) * (1/3) ; 

ft6  =  Fs/2*linspace (0, l,NFFT/2+l) ; 


%  plot  inverse  FFT  x 
Fcx  =  Fx; 

f Fcx  =  NFFT*ifft (Fcx, NFFT) ; 

figure (fig) 
cl  f 

fig  =  fig  +  1 

plot (tt6, x,  ' r ' )  ; 
hold  on 

plot (tt6, f Fcx ( 1 : lenx) ,' g* ') ;  %  !!!  NFFT  ->  true  length 
ylabel ( 'Elevation  [in]  ' )  ; 
xlabel (' Distance  [in]'); 
title (' Munson  Gravel') 

%title (' Perryman  3') 

'size  x, f Fcx ' 

size (x) ,size ( f Fcx) 

plot  (tt6, x-fFcx (1 : lenx) ,  'b ' ) ; 

legend ( ' original ' , ' original (FFT/ iFFT) ' , ' error ' , ' Location ' , ' Northwest ' ) 
print ( ' -dpng ' , ' t6-ter ' ) ; 

figure (fig) 
cl  f 

fig  =  fig  +  1 

plot (tt6, x,  ' r ' )  ; 
hold  on 

plot (tt6 , f Fcx ( 1 : lenx) ,' g* ') ;  %  !!!  NFFT  ->  true  length 
ylabel (' Elevation  [in]'); 
xlabel (' Distance  [in]'); 
title (' Munson  Gravel') 

%title (' Perryman  3') 


'size  x, f Fox ' 
size (x) , size (fFcx) 

legend ( ' original ' ,  ' original (FFT/ iFFT) ' ,  ' Location ' ,  ' Northwest ' ) ; 
axis ( [1.45e4,1.4535e4,1733,1736] ) ; 
print ( ' -dpng ' , ' t6-zm ' ) ; 

%t7  .  m 


%  plot  exp  filter 
H  =  Ffexp; 

%H=conv (Ffexp, Frect) ; 
tmp  =  size (H) ; 
lenH  =  tmp(l); 


%  plot  inverse  x,  mean,  mydet,  detrended 
Fcx  =  Fx.*H; 

%Fcx  =  Fx . * (lenH*H) . *Frect; 

fFcx  =  NFFT*NFFT*if f t (Fcx, NFFT) ;  %  !!!  scale  by  time  length 

figure (fig) 
cl  f 

fig  =  fig  +  1 

plot ( fmean, ' r ' ) 
hold  on 

plot (abs ( fFcx (lenfexp : lenx+ (lenfexp-1 )- (lenfexp-1 ))),' g* ') ;  %  !!!  proper  indexes  for  y  cutoffs  (y 
only  when  exp  fully  in  orig  y) 

'size  fmean, fFcx, lenx, lenfexp, lenN ' 
size (fmean) , size (fFcx) , lenx, lenfexp, lenN 
title (' Munson  Gravel  -  Mean'); 
ylabel (' Elevation  [in]'); 
xlabel (' Distance  [in]'); 

legend ( ' spacial  filter ' , ' frequency  filter (FFT/ . *H/ iFFT) ' , ' Location ' , ' Northwest ' ) ; 
print ( ' -dpng ' , ' t 7 -mean ' ) ; 

figure (fig) 
cl  f 

fig  =  fig  +  1 

plot ( fmean, ' r ' ) 
hold  on 

plot (abs (fFcx (lenfexp : lenx+ (lenfexp-1) - (lenfexp-1) )),' g* ') ;  %  !!!  proper  indexes  for  y  cutoffs  (y 
only  when  exp  fully  in  orig  y) 

'size  fmean, fFcx, lenx, lenfexp, lenN ' 
size (fmean) , size (fFcx) , lenx, lenfexp, lenN 
title (' Munson  Gravel  -  Mean'); 

%title ( 'Perryman  3  -  Mean7 ) ; 
ylabel (' Elevation  [in]'); 

legend ( ' spacial  filter ' , ' frequency  filter (FFT/ .  *H/ iFFT) ' , ' Location ' , ' Northwest ' ) ; 
axis ( [4 .333e4, 4 .3345e4, 1735.25,  1735.45]); 
print ( ' -dpng ' , ' t7-zm ' ) ; 

%t8  .  m 


%  plot  detrended 

%  assumes  lenexp  is  odd,  take  out  mean  from  FFT/.*H/iFFT 

mydet  =  x(  (lenfexp+1 ) /2  :  lenx  -  (lenfexp+1) /2  +  1  )  -  abs (  f Fcx (  lenfexp  :  lenx  +  (lenfexp-1)  - 
(lenfexp-1)  )  ) ; 

mm=mean (mydet) 
mydet  =  mydet  -  mm; 

figure (fig) 
cl  f 

fig  =  fig  +  1 

plot ( fdet ,  ' r ' )  ; 
hold  on 

plot (  mydet, ' g* ' ) 

title (' Munson  Gravel  -  Detrended'); 
ylabel ( 'Elevation  [in]  ' )  ; 


45 


xlabel (' Distance  [in]'); 

legend (' original  -  spacial  f ilter ',' original  -  frequency 
filter ( FFT/ . *H/ iFFT) ' , ' Location ' , ' Northwest ' ) ; 
print ( ' -dpng ' ,  ' t8-det ' ) ; 

figure (fig) 
cl  f 

fig  =  fig  +  1 

plot (fdet, ' r ' ) ; 
hold  on 

plot (  mydet, 'g*') 

title (' Munson  Gravel  -  Detrended'); 
ylabel (' Elevation  [in]'); 
xlabel (' Distance  [in]'); 

legend (' original  -  spacial  f ilter ',' original  -  frequency 
filter (FFT/ . *H/iFFT) ' , ' Location ' , ' Northwest ' ) ; 
axis ( [4.331e4,4.337e4,-0.4,0.8]  )  ; 
print ( ' -dpng ' , ' t8-zm ' ) ; 

%t9  .m 


%  plot  mean  error 

figure (fig) 
cl  f 

fig  =  fig  +  1 

subplot (2 , 1 , 1 ), plot ( fmean  -  abs (fFcx (lenfexp : lenx+ (lenfexp-1) - (lenfexp-1) ) ) ,  ' r ' ) ;  %  !!!  proper 
indexes  for  y  cutoffs  (y  only  when  exp  fully  in  orig  y) 
title (' Munson  Gravel  Spacial/Frequency  Filter  Differences') 
ylabel ( ' Mean ' ) ; 


%  plot  detrended  error 

subplot (2,1,2) , plot ( fdet -mydet ,  'r'); 

%title (' Spacial/Frequency (iFFT)  Munson  Gravel  Differences') 
ylabel ( ' Detrended ' ) ; 
print ( ' -dpng ' , ' p9-dif f ' ) ; 

%t!0  . m 


%  PSD 

%  std(x,l)  =  sqrt (sum (x . A2 ) /length (x) ) 

tmp  =  size (mydet) ; 
lenmydet  =  tmp(l) 

Fmydet  =  f ft (mydet , NFFT) /NFFT; 

Pxx  =  Fmydet con j (Fmydet) ; 

! rm  tlO-d 
diary  tlO-d 
' spacial/freq  restuls' 

'  rms  ' 

' std (mydet, 1 )  ' 

' sqrt (var (mydet, 1 ) )  ' 

' sqrt (  sum  (mydet . A2 ) / length (mydet) )  ' 

' sqrt (  sum (  Pxx  )* length (mydet ) /NFFT) ' 

[rms , std (mydet , 1 ), sqrt (var (mydet , 1 )), sqrt (  sum (mydet . A2 ) /length (mydet )), sqrt (  sum (  Pxx 
) *NFFT/length (mydet) ) ] 

%'Fs/2',  sqrt (  sum (2*Pxx) *Fs  ) 

diary 

stop 

%EXAMPLE :  Spectral  analysis  of  a  complex  signal  plus  noise. 

%  Fs  =  1000;  t  =  0 : 1/Fs : .296; 

%  x  =  exp (li*2*pi*200*t) +randn (size (t) ) ; 

%  h  =  spectrum. periodogram;  %  Create  a  periodogram  spectral  estimator. 

%  psd (h, x, ' Fs ' , Fs) ;  %  Calculates  and  plots  the  two-sided  PSD. 

%  EXAMPLE:  Spectral  analysis  of  a  signal  that  contains  a  200Hz  cosine 
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%  plus  noise. 

Fs  =  1000;  t  =  0 : 1/Fs : .296; 
x  =  cos (2*pi*t*200 ) +randn (size (t) ) ; 
h  =  spectrum. welch; 
psd (h,xf  ' Fs ' , Fs ) ; 


%  Create  a  Welch  spectral  estimator. 
%  Calculate  and  plot  the  PSD. 


APPENDIX  A.5  -  Test  Case  Input/Output  Example 

%tn 

%  %%%!!!  divid  by  NFFT  (assumes  zero  mean,  else  need  add  mean/Pxx (0)  etc.) 
clear  all 
diary  tll-d 

%%%%%!  !  !  !  ! NEED  TO  HAVE  ZERO  MEAN  or  changes  things  for  RMS 
x  =  [1  3  -10  5  8  9  22  2  3  10  10  12  11  14  13  9  17]  '; 
x  =  x-mean (x) ; 
tmp  =  size (x) ; 
lenx  =  tmp ( 1 ) ; 

xl  =  [1  3  -10  5  8  9  22  2  3  10  10  12  11  14  13  9  17]  '; 
xl  =  [ xl ; 1 ; 0 ; - 1 ; 2  0 ; 3  0 ] ; 
xl  =  xl  -mean (xl) ; 

tmp  =  size (xl ) ; 
lenxl  =  tmp(l); 

%%%h 

h  =  [1  2  3  8  3  2  1]'; 
tmp  =  size (h) ; 
lenh  =  tmp ( 1 ) ; 

' nextpow2 ' 
lenx 

nextpow2 (lenx) 

NFFT  =  2Anextpow2 (lenx) 

%%%convolution 
xh=conv (x, h, ' full ' ) ; 

Fx  =  fft (x, NFFT) /NFFT; 

Fh  =  fft  (h, NFFT) /NFFT; 

Fxh  =  Fx . *  Fh ; 

iFxh  =  NFFT*NFFT*ifft (Fxh,NFFT) ; 
iFxh  =  iFxh ( 1 : lenx+lenh-1 ) ; 
iFx  =  NFFT*ifft (Fx,NFFT)  ; 
iFh  =  NFFT*ifft (Fh,NFFT) ; 

'  x,  h ' 

x,  h 

size (x) , size (h) 

'xh' 

xh, iFxh 

'size  xh,iFxh' 

size (xh) , size (iFxh) 

' fft (x, NFFT) /lenx' 

Fx 

'if ft (Fx, NFFT) *lenx ' 
iFx 

' conj  (Fx) ' 
conj  (Fx) 

' var (x) , var  (x, 1 )  ' 
var (x) , var (x, 1 ) 

' std (x) , std (x, 1 ) ' 
std (x) ,  std (x, 1 ) 

' exp (-x) ' 
exp (-x) 

' linspace ( 0 , 1 , NFFT/2+1 ) ' 

NFFT 

NFFT/2+1 

linspace (0, 1, NFFT/2+1) 

'  psd ' 

Fs  =  1/3;  %T=3 

hl=psd (spectrum. welch, 0 . 5* sin (2*pi*Fs* [0:3:1000]), ' Fs ' , Fs) 
fprintf ( ' %f  %e\n ' , [hi . frequencies ' ; hi . data ' ] ) ; 


%  PSD 

Fxl  =  fft (xl, NFFT) /NFFT; 
' var  x, xl ' 
al=var (x, 1 ) 
bl=var (xl , 1 ) 

Fs  =  1/3; 
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Pxx  =  Fx.*conj  (Fx) ; 

Pxl  =  Fxl . *conj  (Fxl ) ; 

'sum  x.A2  Pxx*NFFT;  al  sum  Pxx*NFFT/lenx ' 

[sum(x. A2)  sum (Pxx) *NFFT] 

[al  sum (Pxx) *NFFT/lenx] 

%%%!!!  divid  by  NFFT  (assumes  zero  mean,  else  need  add  mean/Pxx (0)  etc.) 

'sum  Px  Pxl  /NFFT' 

cl=sum (Pxx) *NFFT/ length (x) 

dl=sum (Pxl ) *NFFT/ length (xl ) 

rl  =  (  al*length (x)  )  /  (  cl*length (x)  ); 
r2  =  (  bl*length (xl)  )/  (  dl*length (xl )  ); 

'al*len,  cl*lenA2  ratio,  rep  bl,dl' 

fprintf('%f  %f  %f \n ' ,  [al*length (x) , cl*length (x) ]  ,  rl )  ; 
fprintf('%f  %f  %f \n ' , [bl*length (xl ) , dl* length (xl ) ] , r2 ) ; 

'sqrt(al),  sqrt (cl/length (x) ) ' 
sqrt (al) , sqrt (cl) 

' std (x, 1 ) ' 

std(x,l)  %  sqrt (al*length (x) ) /sqrt (length (x) ) 

' sqrt (var (x, 1) ^length (x) ) ,  std (x, 1) *sqrt (length (x) ) ' 
sqrt (var (x, 1 ) * length (x) ) , std (x, 1 ) *sqrt (length (x) ) 
al* length (x) , std (x, 1 ) A2* (length (x) ) 
diary 


%===tll-d 


ans  = 
nextpow2 


lenx  = 


17 


ans  = 


5 


NFFT  = 


32 


ans  = 
x,  h 


-7.1765 

-5.1765 

-18.1765 

-3.1765 

-0.1765 

0.8235 

13.8235 

-6.1765 

-5.1765 

1.8235 

1.8235 

3.8235 

2.8235 

5.8235 

4.8235 

0.8235 

8.8235 
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h 


1 

2 

3 

8 

3 

2 

1 


ans 

1  7  1 


ans 

7  1 


ans  = 

xh 


xh  = 


-7 

.1765 

-19 

.5294 

-50 

.0588 

-112 

.4706 

-124 

.0000 

-184 

.3529 

-82 

.5294 

-28 

.5294 

5 

.4706 

82 

.4706 

-16 

.5294 

-18 

.5294 

16 

.4706 

26 

.4706 

59 

.4706 

67 

.4706 

89 

.4706 

85 

.  6471 

62 

.0000 

88 

.5294 

32 

.  9412 

18 

.4706 

8 

.8235 

iFxh 


-7  . 

1765 

-19. 

5294 

-50. 

0588 

-112  . 

4706 

-124  . 

0000 

-184  . 

3529 

-82. 

5294 

-28. 

5294 

5. 

4706 

82. 

4706 

-16. 

5294 

-18. 

5294 

16. 

4706 

26. 

4706 

59. 

4706 

67. 

4706 

89. 

4706 

85. 

6471 
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62.0000 

88.5294 

32.9412 

18.4706 

8.8235 


ans  = 

size  xh, iFxh 


ans  = 


23  1 


ans  = 


23  1 


ans 


f 1 t  (x, NFFT) / lenx 


Fx  = 

0.0000 

-1.6054 

-0.4515 

-0.9254 

-0.5466 

0.1333 

0.9118 

0.2412 

-0.0993 

-0.3327 

0.7963 

0.2423 

0.1605 

-1.1850 

-0.4037 

-0.5682 

0.0882 

-0.5682 

-0.4037 

-1.1850 

0.1605 

0.2423 

0.7963 

-0.3327 

-0.0993 

0.2412 

0.9118 

0.1333 

-0.5466 

-0.9254 

-0.4515 

-1.6054 


ans 


0 . 1 955i 
+  0 . 8507i 
+  0 . 2054i 
+  1 . 1821i 
+  1 . 107  9i 
+  0 . 4100i 

-  0 . 3004i 

-  0 . 2500i 
+  0 . 4  903i 
+  0 . 1114i 

-  0 . 4 7 7 9i 

-  1 . 0054i 

-  0 . 4 994i 
+  0.177H 
+  0 . 1820i 

-  0 . 1820i 

-  0.177H 
+  0 . 4  994i 
+  1 . 0054i 
+  0 . 4  7  7  9i 

-  0 . 1114i 

-  0 . 4 903i 
+  0 . 2500i 
+  0 . 3004i 

-  0 . 4100i 

-  1 . 107  9i 

-  1.182H 

-  0 . 2054i 

-  0 . 8507i 
+  0 . 1 955i 


if ft (Fx,NFFT) *lenx 


iFx  = 


-7.1765 

-5.1765 

-18.1765 

-3.1765 
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-0.1765 

0.8235 

13.8235 

-6.1765 

-5.1765 

1.8235 

1.8235 

3.8235 

2.8235 

5.8235 

4.8235 
0.8235 

8.8235 
0.0000 

0 

0.0000 

0.0000 

0 

-0.0000 

0.0000 

0 

-0.0000 

-0.0000 

0.0000 

-0.0000 

-0.0000 

0.0000 

0 


ans  = 
conj  (Fx) 


ans 


0 . 0  0  0  0 
-1.6054  + 
-0.4515  - 
-0.9254  - 
-0.5466  - 
0.1333  - 
0.9118  - 
0.2412  + 
-0.0993  + 
-0.3327  - 
0.7963  - 
0.2423  + 
0.1605  + 
-1.1850  + 
-0.4037  - 
-0.5682  - 
0.0882 
-0.5682  + 
-0.4037  + 
-1.1850  - 
0.1605  - 
0.2423  - 
0.7963  + 
-0.3327  + 
-0.0993  - 
0.2412  - 
0.9118  + 
0.1333  + 
-0.5466  + 
-0.9254  + 
-0.4515  + 
-1.6054  - 


ans  = 


0 . 1 955i 
0 . 8507i 
0 . 2054i 
1.182H 
1 . 107  9i 
0 . 4100i 
0 . 3004i 
0 . 2500i 
0 . 4  903i 
0 . 1114i 
0 . 4  7  7  9i 
1 . 0054i 
0 . 4  994i 
Q.177H 
0 . 1820i 

0 . 1820i 
0.177H 
0 . 4  994i 
1 . 0054i 
0 . 4  7  7  9i 
0 . 1114i 
0 . 4  903i 
0 . 2500i 
0 . 3004i 
0 . 4100i 
1 . 107  9i 
1.18211 
0 . 2054i 
0 . 8507i 
0 . 1 955i 


var (x) , var (x,  1 ) 
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ans 


52.5294 


ans  = 


49.4394 


ans  = 

std (x) , std (x,  1 ) 


ans  = 


7.2477 


ans  = 


7.0313 


ans  = 
exp (-x) 


ans  = 


1.0e+07  * 

0.0001 

0.0000 

7.8332 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 


ans  = 

linspace (0,1, NFFT/ 2+1 ) 


NFFT  = 


32 


ans  = 


17 


ans  = 

Columns  1  through  12 


53 


0 

5625  0. 

0.0625  0.1250 

6250  0.6875 

0.1875 

0.2500 

0.3125 

0.3750 

0.4375 

0.5000 

Columns  13 

through  17 

0.7500 

0.8125  0.8750 

0.9375 

1.0000 

ans  = 
psd 


hi  = 


Name : 
Data : 
SpectrumType : 
Normal izedFrequency : 

Fs : 

Frequencies : 
ConfLevel : 
Conf Interval : 


'Power  Spectral  Density' 
[129x1  double] 

' Onesided ' 

false 

0.3333 

[129x1  double] 

'Not  Specified' 

[] 


0.000000 

0.001302 

0.002604 

0.003906 

0.005208 

0.006510 

0.007812 

0.009115 

0.010417 

0.011719 

0.013021 

0.014323 

0.015625 

0.016927 

0.018229 

0.019531 

0.020833 

0.022135 

0.023438 

0.024740 

0.026042 

0.027344 

0.028646 

0.029948 

0.031250 

0.032552 

0.033854 

0.035156 

0.036458 

0.037760 

0.039062 

0.040365 

0.041667 

0.042969 

0.044271 

0.045573 

0.046875 

0.048177 

0.049479 

0.050781 

0.052083 

0.053385 

0.054688 

0.055990 

0.057292 

0.058594 

0.059896 

0.061198 

0.062500 

0.063802 


5 . 22  6231e-25 
9 . 483823e-25 
7 . 050443e-25 
4 . 228587e-25 
1 . 984093e-25 
6 . 87  6082e-2  6 
1 . 569868e-26 
1.860584e-27 
2 . 391157e-28 
2 . 802  668e-28 
2.28592  6e-28 
2 . 048480e-28 
2 . 34  6595e-28 
3 . 117  960e-28 
3. 945803e-28 
4.27  4355e-28 
4 . 732021e-28 
7 . 948712e-28 
2 . 104163e-27 
5 . 810612e-27 
1 . 334  4  42e-2  6 

2 . 4  4 122  9e-2  6 
3 . 586796e-26 

4 . 2  96287e-2  6 
4 . 240851e-2  6 

3 . 4  64528e-2  6 
2 . 338243e-2  6 

1 . 2  9907  6e-2  6 
6 . 011988e-27 
2 . 513154e-2  7 
1 . 17  6740e-27 

7 . 2  97  665e-28 
5 . 569431e-28 
4 . 908393e-28 
4 . 682606e-28 
4 . 53232  9e-28 
4 . 708755e-28 
5 . 515843e-28 
6 . 72  6700e-28 
8 . 231078e-28 
1 . 109477e-27 
1 . 771123e-27 
3 . 1054  92e-27 
5 . 3352  7  7e-2  7 
8 . 395125e-27 
1 . 172202e-2  6 
1 . 430514e-2  6 
1 . 513542e-2  6 
1 . 380  611e-2  6 
1 . 07  9887e-2  6 
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0.065104 

0.066406 

0.067708 

0.069010 

0.070312 

0.071615 

0.072917 

0.074219 

0.075521 

0.076823 

0.078125 

0.079427 

0.080729 

0.082031 

0.083333 

0.084635 

0.085938 

0.087240 

0.088542 

0.089844 

0.091146 

0.092448 

0.093750 

0.095052 

0.096354 

0.097656 

0.098958 

0.100260 

0.101562 

0.102865 

0.104167 

0.105469 

0.106771 

0.108073 

0.109375 

0.110677 

0.111979 

0.113281 

0.114583 

0.115885 

0.117188 

0.118490 

0.119792 

0.121094 

0.122396 

0.123698 

0.125000 

0.126302 

0.127604 

0.128906 

0.130208 

0.131510 

0.132812 

0.134115 

0.135417 

0.136719 

0.138021 

0.139323 

0.140625 

0.141927 

0.143229 

0.144531 

0.145833 

0.147135 

0.148438 

0.149740 

0.151042 

0.152344 

0.153646 

0.154948 

0.156250 

0.157552 

0.158854 

0.160156 


7 . 214590e-27 
4 . 147154e-27 
2 . 151817e-27 
1 . 147077e-27 

7 . 2  4  03  67e-2  8 

5 . 2  96537e-28 
4 . 161614e-28 
3.54771 9e-2  8 
3 . 287241e-28 
3 . 259712e-28 
3 . 57  667  9e-28 
4 . 289426e-28 
5 . 180134e-28 
6 . 2537  69e-28 
8 . 344877e-28 
1 . 289231e-27 
2 . 090314e-27 
3 . 191044e-27 
4 . 360914e-27 
5 . 248827e-27 
5 . 538457e-27 
5 . 117029e-27 
4 . 136936e-27 
2 . 923833e-27 
1 . 809939e-27 
1 . 004219e-27 
5 . 478430e-28 

3 . 4  934  97e-28 
2 . 738669e-28 
2 . 308483e-28 
1 . 966835e-28 
1.74779  6e-2  8 
1 . 634736e-28 
1 . 669708e-28 
2 . 0432  68e-28 
2 . 827591e-28 
3 . 816106e-28 
4 . 964  963e-28 
7 . 017416e-28 
1 . 141482e-27 
1 . 926689e-27 
3 . 0114  99e-27 
4 . 14  9922e-27 
4 . 983957e-27 
5 . 222  663e-27 

4 . 7  94  617e-27 
3 . 867756e-27 
2 . 739608e-27 
1 . 696213e-27 
9 . 218688e-28 
4 . 6912  63e-28 
2 . 750187e-28 
2 . 2107  99e-28 
2 . 098549e-28 
2 . 009576e-28 
1 . 903970e-28 
1 . 822989e-28 
1 . 894325e-28 
2 . 360866e-28 
3 . 378028e-28 
4 . 905343e-28 
7 . 012244e-28 
1 . 019448e-27 
1 . 505093e-27 
2 . 142128e-27 

2 . 7  857  4  8e-2  7 
3 . 212215e-27 
3 . 24  9653e-27 
2 . 890991e-27 

2 . 2  937  55e-2  7 
1 . 664718e-27 
1 . 135004e-27 
7 . 330027e-28 

4 . 4  4  94  90e-2  8 


0.161458 

0.162760 

0.164062 

0.165365 

0.166667 


2 . 637  662e-28 
1 . 785309e-28 
1 . 55157  9e-28 
1 . 5392  91e-28 
7 . 735081e-2  9 


ans  = 


ans  = 
var  x,xl 

al  = 


49.4394 


bl  = 


75.2417 


ans  = 

sum  x.A2  Pxx*NFFT;  al  sum  Pxx*NFFT/lenx 


ans  = 

840.4706  840.4706 


ans  = 


49.4394  49.4394 


ans  = 

sum  Px  Pxl  /NFFT 


cl  = 


49.4394 


dl  = 


75.2417 


ans  = 

al*len,  cl*lenA2  ratio,  rep  bl,dl 

840.470588  840.470588  1.000000 
1655.318182  1655.318182  1.000000 

ans  = 

sqrt(al),  sqrt (cl/length (x) ) 

ans  = 

7.0313 


ans 


56 


7.0313 


ans  = 
std (x,  1 ) 


ans 


'/ .  0  3  1  3 


ans  = 

sqrt (var (x, 1) ^length (x) )  ,  std (x, 1 

ans  = 

28.9909 

ans 

3  8.9  9  3  9 

ans  = 

848.4706 

ans  = 

849.4706 


)  *sqrt (length (x) ) 


57 


